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Abstract 

Compressible  stability  of  growing  boundary  layers  is  studied  by  numerically 
solving  the  partial  differential  equations  under  a  parabolizing  approximation.  The 
resulting  parabolized  stability  equations  (PSE)  account  for  non-parallel  as  well  as 
nonlinear  effects.  Evolution  of  disturbances  in  compressible  flat-plate  boundary 
layers  are  studied  for  freestream  Mach  numbers  ranging  from  0  to  4.5.  Results  indi¬ 
cate  that  the  effect  of  boundary-layer  growth  is  important  for  linear  disturbances. 
Nonlinear  calculations  are  performed  for  various  Mach  numbers.  Two-dimensional 
nonlinear  results  using  the  PSE  approach  agree  very  well  with  those  from  direct  nu¬ 
merical  simulations  using  the  full  Navier-Stokes  equations  while  the  required  com¬ 
putational  time  is  less  by  an  order  of  magnitude.  Spatial  simulations  using  PSE 
have  been  carried  out  for  both  the  fundamental  and  subharmonic  type  breedcdown 
for  a  Mach  1.6  boundary  layer.  The  promising  results  obtained  in  this  study  show 
that  the  PSE  method  is  a  powerful  tool  for  studying  boundary-layer  instabilities 
and  for  predicting  transition  over  a  wide  range  of  Mach  numbers. 
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I.  Introduction 


The  subject  of  compressible  boxmdary-layer  stability  has  attracted  a  great  deal 
of  interest  in  the  past  few  years  due  to  its  importance  in  understanding  the  on¬ 
set  of  transition  in  high-speed  flows  and  providing  some  theoretical  background 
for  laminar  flow  control  (LFC)  techniques  (Malik,  1990a).  Most  investigations  of 
compressible  linear  stability  (e.g.,  Mack,  1969,  1984)  have  employed  what  is  known 
as  the  “quasi-parallel”  approach  whereby  the  growth  of  the  boundary  layer  is  ig¬ 
nored  and  the  linearized  Navier-Stokes  equations  are  reduced  to  ordinary  differential 
equations  (ODE)  by  assuming  a  wave-like  disturbance  of  the  form 

(1) 

where  x,  y,  and  z  are  the  streamwise,  wall-normal,  and  spanwise  coordinates,  respec¬ 
tively;  o  and  yS  are  the  corresponding  wave  numbers,  a;  is  the  disturbance  frequency 
and  ’J'  represents  the  disturbance  shape  function.  The  linear  ODE’s  along  with  the 
homogeneous  boundary  conditions  constitute  an  eigenvalue  problem  of  the  form 

a  =  a(w,  $)  (2) 

which  can  be  solved  by  standard  eigenvalue  techniques.  The  imaginary  part  of 
a  gives  the  disturbance  growth  rate  and  a  small  disturbance  is  expected  to  grow 
provided  Oj  <  0.  For  a  given  flow,  this  eigenvalue  approach  can  be  applied  “locally” 
at  various  locations  along  the  body  in  order  to  obtsdn  an  idea  about  overall  growth 
of  disturbances  and  to  correlate  with  transition  location  using  empirical  methods 
such  as  the  method. 

The  effect  of  non-parallel  flow  on  boundary-layer  instability  has  been  studied 
by  Gaster  (1974),  Saric  and  Nayfeh  (1975),  Gaponov  (1981),  and  El-Hady  (1991). 
In  the  multiple-scales  method  used  by  the  latter  three  authors,  the  disturbances 
are  decomposed  into  a  slowly- varying  shape  function  and  a  rapidly-oscillating  wave 
part.  Both  parts  are  represented  as  functions  of  a  fast-scale  variable  (x)  and  a  slow- 
scale  variable  (x  =  ex,  with  e  =  1/i?).  With  these  assumptions  the  governing  PDE’s 
are  reduced  to  a  set  of  ODE’s  by  neglecting  terms  of  order  equal  to  or  higher  than 
c^.  In  conjunction  with  the  solvability  condition,  the  analysis  yields  non-parallel 
corrections  to  the  eigenvalues  computed  by  the  quasi-parallel  theory.  Just  like  the 
traditional  linear  theory,  the  multiple- scales  approach  can  only  be  applied  locally 
for  a  given  problem. 

Apart  from  the  “local”  methods  described  above,  the  evolution  of  disturbances 
in  a  given  flowfield  may  also  be  computed  numerically  by  solving  the  governing 
partial  differential  equations  (PDE’s)  without  resorting  to  the  eigenvalue  approach. 
The  effect  of  boundary-layer  growth  and  other  history  effects  associated  with  ini¬ 
tial  conditions  and  variation  in  wall  temperature,  for  instance,  can  be  properly 
accounted  for.  This  was  done  for  the  Gortler  vortex  problem  by  Hall  (1983)  and 
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Spall  and  Malik  (1989).  Denier  et  al.  (1991)  solved  the  “receptivity”  problem  to 
provide  the  inflow  conditions  for  the  PDE’s  and  were  able  to  show  how  Gortler 
vortex  structure  develops  from  a  discrete  roughness  site. 

The  governing  PDE’s  for  the  Gortler  problem  are  parabolic  and  thus  the  solu¬ 
tion  can  be  obtained  by  direct  marching  provided  the  initial  conditions  are  known. 
However,  the  governing  equations  for  Tollmien-Schhchting  (TS)  and  inviscid  type 
disturbances  are  elliptic  and  their  solution  cannot  be  obtained  by  simple  march¬ 
ing  methods.  In  addition,  the  nximerical  solution  of  these  PDE’s  requires  proper 
outflow  boimdary  conditions  which  is  a  nontrivial  task.  However,  we  note  that  for 
boundary-layer  type  flows  which  are  of  interest  here,  the  equation  set  is  only  weakly 
elliptic  along  the  dominant  flow  direction.  Therefore,  with  appropriate  simplifica¬ 
tions,  one  could  “parabolize”  these  stability  equations  and  avoid  the  difficulties 
associated  with  the  downstream  boundary  conditions. 

From  a  physical  view  point,  the  streamwise  ellipticity  arises  from  the  upstream 
propagation  of  acoustic  waves  and  the  streamwise  viscous  diffusion.  To  render 
the  stability  equations  parabolic,  one  must  devise  a  way  to  suppress,  but  without 
compromising  the  essential  physics,  this  upstream  propagation.  One  way  to  derive 
the  parabolized  stability  equations  (PSE)  is  to  borrow  ideas  &om  the  multiple-scales 
approach  and  decompose  the  disturbance  into  a  rapidly- varying  wave-like  part  and  a 
slowly- varying  shape  function.  The  ellipticity  is  retained  for  the  wave  part  while  the 
peurabolization  is  applied  to  the  shape  fimction.  The  resulting  PSE  can  be  solved  by 
marching  along  the  streamwise  direction.  The  technique  can  be  used  to  study  both 
the  linear  and  nonlinear  evolution  of  convective  distiirbances  in  growing  boimdary 
layers.  Global  or  absolute  instabilities  can  not  be  studied  by  this  approach.  This 
parabolizing  procedure  has  been  used  recently  by  Bertolotti  et  al.(1992)  for  Blasius 
flow. 

The  objective  of  this  research  is  to  study  compressible  boimdary-layer  stability 
and  transition.  We  employ  parabolized  stability  equations  for  linear  and  nonlin¬ 
ear  development  of  disturbances  in  a  compressible  boimdary  layer.  The  nonlinear 
calculations  are  carried  all  the  way  to  the  transition  stage  for  supersonic  flows.  In 
section  II,  we  formulate  the  problem  while  the  numerical  procedure  used  to  solve 
the  governing  equations  are  given  in  section  III.  The  results  and  conclusions  are 
presented  in  section  FV  and  V,  respectively. 


II.  Problem  Formulation 

The  evolution  of  disturbances  in  compressible  boundary  layers  is  governed  by 
the  compressible  Navier-Stokes  equations 

^  +  V.(pK)  =  0 
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(3) 


p[^  +  (y  ■  v)v]  =  -  vp  -1-  v(A(v  •  F)]  +  V  •  [p(  Vi?  +  vi?^)] 
pcp[^  +  (V  •  V)T]  =  V .  (kVT)  +  ^  +  (i?  .  V)p  +  $ 

where  V  is  the  velocity  vector,  p  the  density,  p  the  pressure,  T  the  temperature,  Cp 
the  specific  heat,  k  the  thermal  conductivity,  p  the  first  coefficient  of  viscosity,  and 
A  the  second  coefficient  of  viscosity.  The  viscous  dissipation  function  is  given  as 

$  =  A(V  •  V)^  +  ^[VV  +  Vi?^]2. 

The  equation  of  state  is  given  by  the  perfect  gas  relation 

p  =  p3?T 


and  the  steady  state  solution  of  the  basic  flow  can  be  derived  by  invoking  the 
boundary-layer  assumption. 

In  this  research,  we  formulate  the  compressible  stability  problem  in  Cartesian 
coordinates  for  the  flat-plate  geometry,  although  the  theory  itself  can  be  easily 
extended  to  axisymmetric  bodies  amd  infinite  swept-wing  flows.  The  Cartesian 
coordinates  are  denoted  by  r,  y,  and  z  to  represent  the  streamwise,  wall-normal, 
and  spanwise  directions,  respectively.  All  the  lengths  are  scaled  by  a  reference  length 
I,  velocity  by  u*,  density  by  p*,  pressure  by  Pe«e»  time  by  //u*,  and  other  variables 
by  the  corresponding  boimdary-layer  edge  values.  The  basic  flow  is  perturbed  by 
fluctuations  in  the  flow,  i.e.  the  total  field  can  be  decomposed  into  a  mean  value 
(boundary-layer  solution)  and  a  perturbation  quantity 

u  =  u-\-u\  V  =  V  +  v',  w  =  w  +  w' 

p  =  p  +  p',  p  =  p-|-p',  T  =  f  +  T'  (4) 

p  —  p,  p' ,  A  =  A  A^,  k  =  k  k\ 


Substituting  Eq.(4)  into  the  Navier-Stokes  equations  given  by  Eq.  (3)  md  sub¬ 
tracting  from  the  governing  equations  corresponding  to  the  steady  mean  flow,  and 
using  the  equation  of  state,  we  obtain  the  governing  equations  for  the  disturbances 
as 


I  yxy 


dxdy 


yy 


+  V, 


ay 

dxdz 


+  v: 


yz 


ay 

dydz 


+  V, 


dU 

dz^ 


(5) 


where  contains  the  disturbance  vector  and  is  defined  as 


(f>  =  {p\u\v',w',T')'^. 
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Matrices  F,  A,  B,  C,  D,  V**,  V*^,  Vyy,  Fu,  Vyt,  and  F„  are  Jacobians  of  the 
corresponding  total  flux  vectors  and  are  composed  of  a  linear  part  with  only  mean 
flow  quantities  (denoted  by  superscripts  1)  and  a  nonlinear  part  which  contains 
perturbation  quantities  (denoted  by  superscripts  n):  F  =  F*  +  F”,  A  =  A*  +  A”, 
etc.  We  note  here  that  matrices  F,  A,  C,  D  have  contributions  from  both  inviscid 
and  viscous  terms,  and  thus  contain  terms  of  order  one  and  of  order  1/Rq  (Bo  is 
the  reference  Reynolds  number  Ro  =  while  matrices  Fn,  Fu,  V^y,  Vyy, 

^xzi  Vyzi  and  Vzx  are  solely  due  to  viscous  diffusion  and  are  of  order  I/Rq. 

To  facilitate  our  discussion  on  the  relation  between  linear  and  nonlinear  dis¬ 
turbances,  we  rearrange  Eq.  (5)  in  the  following  form 


at  ox  ay  dz 


”5x2  ''^Vdxdy 
’^dxdz  '•^dydz 


»»  5y2 

_  F'  ^  -  F" 
”5z2 


(6) 

where  the  left  hand  side  contains  only  linear  operators  operating  on  the  disturbance 
vector  and  the  right-hand-side  forcing  vector  F”  is  due  to  nonlinear  interaction  and 
includes  all  nonlinear  terms  zissociated  with  the  disturbances.  The  right  hand  side 
is  given  as 

F”  =-F”^-A"^-B”' 
at  ox 


5x2 


5y  dz 
^dxdy  *'5'5y2 


(7) 


+  K.” 


d'^4> 


d^<f> 


+  F” 

”5x5x  ^^dydz 


-L  - 

^  ”5^2 


In  the  incompressible  limit,  F"  contains  quadratic  nonlinearities;  while,  for  com¬ 
pressible  flows,  cubic  and  higher-order  nonlinearities  are  present.  For  small  distur¬ 
bances,  F”  can  be  neglected  and  thus  Eq.  (6)  reduces  to  the  linearized  Navier-Stokes 
equations 

dt 


+  X'^  +  B'^+C'^  +  Z)V 

5x  dy  dz 


-VI 


-  F' 


d^<t> 

-  f'  ^  ^ 

5x2  ■ 

’’didy 

5^</> 

V' 

dydz 

”  5x2 

-F 


yy 


dy^ 


dxdz 


(8) 


=  0. 


The  governing  PDE’s  of  the  disturbances.  Equation  (6),  is  hyperbolic  in  time 
for  the  convection  terms  (inviscid  part).  When  we  consider  only  the  spatied  deriva¬ 
tives,  Equation  (6)  is  elliptic  in  the  streamwise  direction  due  to  two  reasons.  First, 
the  streamwise  viscous  term  Fi*  zdlows  any  disturbances  to  be  diffused  upstream. 
Second,  and  more  importantly,  the  convection  term  in  the  streamwise  direction 


makes  the  upstream  propagation  of  acoustic  waves  possible.  The  latter  can  be  bet¬ 
ter  understood  by  considering  the  linearized  version  of  the  inviscid  equations  and 
using  the  method  of  characteristics  (MOC)  theory.  Since  the  inviscid  part  of  Eq. 
(8)  is  hyperbolic  in  time,  the  corresponding  slopes  of  characteristic  lines  in  the  x  —  t 
plane  (which  determine  the  direction  of  propagation)  can  be  found  by  solving  the 
following  eigenvalue  equation 


\A‘  -  AcT'I  =  0. 


Negative  eigenvalues  imply  the  wave  is  propagating  from  downstream  to  upstream 
and  vice  versa.  The  eigenvalues  of  the  above  equation  are 

Ac  =  u,  ti,  u,  u  -j-  c,  u  —  c 

where  c  is  the  speed  of  sound.  For  boundary-layer  flows  of  interest  in  this  study,  the 
first  four  eigenvalues  are  always  positive,  while  the  last  eigenvalue  (ti  —  c)  can  be 
either  negative  or  positive  depending  upon  the  local  Maeh  number  (M*  =  u/c).  For 
subsonic  flows,  this  quantity  is  negative  throughout  the  whole  flowfield,  therefore, 
the  equation  set  (8),  and  thus  (6),  is  elliptic.  For  supersonic  flows,  the  ellipticity 
only  arises  inside  the  subsonic  layer  adjacent  to  the  wall. 

Based  upon  the  above  discussion,  one  way  to  “parabolize”  the  PDF’s  given  by 
Eq.  (6)  and  make  the  marching  solution  feasible  is  to  neglect  the  viscous  diffusion 
terms  along  the  streamwise  direction  and  prohibit  the  upstream  wave  propagation 
either  by  dropping  the  left-running  ch^acteristics  (associated  with  the  eigenvalue 
u— c)  (Chang  and  Merkle,  1989)  or  suppressing  some  part  of  the  streamwise  pressure 
gradient,  as  it  is  done  in  the  Parabolized  Navier-Stokes  (PNS)  approach  (Vigneron 
et  al.,  1978).  For  the  stability  equations,  the  upstream  wave  propagation  can  be 
suppressed  by  either  dropping  the  characteristic  equation  associated  with  the  eigen¬ 
value  u  —  c  or  multiplying  the  streamwise  pressure  disturbance  gradient  dp'  jdx  by 
a  parameter  ft  given  by 


f2  = 


l-h(7-l)Mr 


Mr  <  1 


(9) 


=  1, 


Mr  >  1 


where  7  is  the  ratio  of  specific  heats.  These  parabolizing  procedures  are  quite 
effective  for  the  PNS  approach  and  yield  solutions  which  compare  favorably  with 
those  obtained  by  the  full  Navier-Stokes  equations  provided  a  large  portion  of  the 
flow  is  supersonic  and  only  steady  state  solutions  are  of  interest  (Vigneron  et  al., 
1978;  Rubin,  1981).  The  advantage,  of  course,  is  the  significant  reduction  in  the 
computational  cost  due  to  the  marching  solution. 

For  compressible  stability  problems,  the  disturbances  are  essentially  unsteady 
waves  propagating  across  the  whole  boundary  layer  and  the  amplitudes  of  these 
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waves  reach  their  maxima  near  the  critical  layer  located  between  the  wall  and  the 
boundary-layer  edge.  These  instability  waves  undergo  a  “fast  oscillation”  (phase 
change)  as  they  evolve  along  the  flow  direction.  Direct  application  of  the  parab¬ 
olizing  procedure  used  in  the  PNS  approach  for  mean  flow  computations  to  our 
governing  stability  equations  would  not  capture  the  flow  physics  due  to  the  sup¬ 
pression  of  the  wave  propagation  along  the  left-running  characteristics.  Therefore, 
an  alternative  procedure  must  be  devised. 

Linear  PSE 

As  mentioned  previously,  one  way  to  “parabolize”  the  governing  PDE’s  is  to 
first  decompose  the  disturbances  into  a  fast-oscillatory  wave  part  and  a  slowly- 
varying  shape  function.  We  keep  the  ellipticity  for  the  wave  part  while  parabolizing 
the  governing  equation  for  the  shape  function.  Following  the  lead  of  the  non-parallel 
linear  stability  theory,  we  assume  that  the  distmrbance  vector  ^  for  an  instability 
wave  with  a  frequency  u  and  a  spanwise  wave  number  ^  (assume  the  wave  is  periodic 
in  both  the  temporal  and  spanwise  directions)  can  be  expressed  as 

(10) 

where  x  is  the  fast-scale  variable,  a(x)  is  the  corresponding  streamwise  wave  number 
and  ^  is  the  “shape  function”  vector  given  by 


=  (p,u,w,t&,  T)^. 


(11) 


As  compared  to  Eq.  (1),  the  shape  function  ^  is  now  a  function  of  both  x  and  y 
due  to  the  growth  of  the  boundary  layer  and  the  wave  number  a  is  a  fimction  of  x 
to  accoimt  for  the  growing  boundary  layer.  For  simplicity,  we  now  restrict  ourselves 
to  the  linear  case,  i.e.,  only  a  single  disturbance  mode  (u;,/3)  is  considered  and  the 
nonlinear  effect  will  be  included  later  on.  Substituting  Eq.  (10)  into  the  linear 
stability  equation  (8),  we  have  the  following  equation  for  the  shape  function 


(12) 


where  the  vectors  D,  A  and  B  are  defined  by 


D  =  -iuT^  +  D*  +  iaA*  -I- 

A  =  A‘-2iaVi,-i0Vi, 
B  =  B'  -  iaVj,  - 
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In  the  quasi-paxallel  linear  theory  where  “normal-mode”  analysis  is  employed, 
the  shape  function  4'  is  assumed  to  be  a  function  of  y  only  (d^’/di  =  0);  therefore, 
Equation  (12)  reduces  to  the  following  system  of  ODE’s 


lo®  =  0 


(13) 


where  the  operator  Lq  is  given  by 


L„=D^B±-Vi,^ 

and  the  elements. of  matrices  D,  B  and  are  evaluated  by  assuming  parallel 
mean  flows  (v  =  0  and  da/dx  =  0).  The  above  ODE’s  in  conjunction  with  homoge¬ 
neous  boundary  conditions  then  constitute  an  eigenvalue  problem  described  by  the 
dispersion  relation  given  in  Eq.  (2). 

Unlike  the  normcil-mode  analysis  described  above,  the  decomposition  (10)  ex¬ 
hibits  some  extent  of  non-uniqueness  between  the  distribution  of  the  wave  part  and 
the  shape  function  part.  In  the  PSE  approach,  we  choose  a  complex  wave  num¬ 
ber  a  and  construct  a  decomposition  such  that  the  change  of  shape  function  'I' 
along  the  streamwise  direction  x  is  of  order  1/Ro  and  the  second  derivative  of  ^ 
{d^’^/dx'^)  is  negligible.  With  this  assumption  and  after  neglecting  all  terms  of 
0(l/jRo),  Equation  (12)  reduces  to 


D<i!  +  A - h  B —  =  V‘  - 

dx  dy  dy^ 


(14) 


Equation  (14)  describes  the  evolution  of  the  shape  function  and  is  “nearly” 
parabolic  in  the  sense  that  second  derivatives  in  x  are  absent  and  the  elliptic  effect 
associated  with  the  wave  part  is  absorbed  in  matrices  D,  A  and  B.  For  instance, 
the  disturbance  pressure  gradient  dp' /dx,  which  is  responsible  for  the  upstream 
influence,  can  be  written  as 

dp'  ..  -  ,  »(/*  oi{x)dx+0z-wt) 

—  =  (.ap+^)eJ: 


The  contribution  of  the  wave  part  (iap)  is  absorbed  in  the  source  term  D'^  and  does 
not  contribute  to  the  upstream  influence  of  the  governing  equations  of  the  shape 
functions,  Eq.  (14).  However,  the  pressure  gradient  shape  function  dp/dx  associ¬ 
ated  with  the  left-running  characteristic  (for  subsonic  flows  only)  is  still  present  in 
the  X  derivative  term.  The  existence  of  this  term  allows  upstream  influence  in  Eq. 
(14).  For  stationary  Gortler  vortex  problem;  a  =  0  and  dp/dx  drops  out,  Eq.  (14) 
reduces  to  the  parabolic  equations  solved  by  Spall  and  Malik  (1989). 

For  supersonic  boundary  layers,  a  large  portion  portion  of  the  flow  possesses 
only  downstream  characteristics,  our  numerical  results  have  shown  that  with  a 
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properly  chosen  value  of  a  (see  discussion  below)  most  of  the  upstream  influence 
is  accounted  for  in  iap  and  the  elliptic  effect  associated  with  the  pressure  gradient 
shape  function,  dp/dx,  is  negligible.  To  make  Eq.  (14)  truly  parabolic  and  enable 
a  stable  marching  procedure  for  subsonic  flows,  we  multiply  dp/dx  by  a  constant 
Q  defined  in  (9)  (in  the  incompressible  limit,  this  is  equivalent  to  setting  dp/dx  to 
zero).  While,  this  is  formally  true  only  in  special  cases  (e.g.  Gortler  vortex  problem), 
the  approximation  yields  solutions  which  compare  very  well  with  accurate  results 
from  full  Navier-Stokes  equations  (Joslin  et  eJ.,  1992).  This  is  because  most  of  the 
ellipticity  is  captured  in  the  iap  term. 

For  incompressible  flows,  one  can  use  the  vorticity-streamfunction  formulation 
for  two-dimensional  flows  or  use  other  formulations  derived  by  eliminating  the  pres¬ 
sure  from  the  momentum  equations,  as  is  done  by  Bertolotti  et  al.  (1992).  In  these 
approaches,  neglecting  second  and  higher  strezimwise  derivatives  of  the  dependent 
variables  inherently  suppresses  some  part  of  the  streamwise  pressure  gradient,  and 
consequently  prohibits  the  upstream  propagation  of  information. 

We  now  describe  the  strategy  to  update  the  streamwise  wave  number  in  or¬ 
der  to  make  the  marching  scheme  well-posed.  The  evolution  of  shape  functions  is 
monitored  during  the  process  of  marching  and  a  is  updated  by  local  iterations  at  a 
given  X  according  to  the  change  in  The  updating  procedure  is  described  herein. 
At  a  given  location  xj ,  we  assume  that  the  streamwise  wave  number  is  given  by  oj 
and  the  total  disturbance  in  the  vicinity  of  xi  can  be  expressed  as 


(f>(x,y,z,t)  =  ’J'(x,  y)e'^'^*>  atdx+0z  wt) 


(15) 


The  change  of  the  shape  function  ^  can  be  approximated  by  the  following  Taylor 
series  expansion  truncated  to  the  first  order 

T/  ^  T 

^(x,y)  =  -I-  -  ^i)  + 

where  is  the  shape  function  at  x  =  xi.  To  an  accuracy  of  0{x  —  xi),  the  above 
equation  can  be  further  expressed  as 


4'(x,y)  = 


(16) 


Substituting  (16)  into  (15),  we  have  the  “effective”  wave  number  in  the  vicinity  of 
xi  given  by 

.  1 

a  =  ai  —  t- - — 

dx 

The  real  part  of  this  effective  wave  number  represents  the  phase  change  of  the 
disturbance  while  the  imaginary  part  depicts  the  growth  rate,  both  corresponding 
to  the  quantity  chosen.  A  disturbance  (^i)  is  unstable  if  the  imaginary  part 
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is  less  than  zero.  The  updating  procedure  of  o  is  repeated  by  using  (17)  until  the 
change  in  o  is  smaller  than  a  prescribed  tolerance  (typically  10“^^). 

Since  the  shape  function  vector  ^  is  a  function  of  y  and  contains  five  depen¬ 
dent  variables  (p,u,  etc.),  the  updating  procedure  above  is  equivalent  to  choosing 
a  normalization  of  the  disturbance  vector  such  that  d^i/dx  is  zero  at  a  particular 
y  location.  Accordingly,  the  value  of  a  computed  by  (17)  will  depend  on  the  y 
coordinate  and  the  selected  dependent  variable  i .  In  this  study,  we  have  used  the 
shape  function  u  (or  t  for  compressible  flows)  at  various  y  locations  or  the  energy 
integral  {E  =  p(u^  -f-  -|-  w^)dy),  which  is  independent  of  the  y  coordinate, 

to  update  the  wave  number  a  and  the  resulting  non-paurallel  growth  rate  (which 
also  depends  on  the  dependent  variable  and  y  coordinate  chosen  to  measure  the 
growth  rate)  appears  to  be  very  weakly  dependent  upon  the  normalization  chosen 
(see  discussion  in  section  IV). 

The  solution  of  (14)  requires  proper  boundary  conditions  in  the  wall- normal 
direction.  We  apply  the  homogeneous  Dirichlet  conditions 

u  =  t;  =  it;  =  r  =  0,  y  =  0  (18) 

at  the  wall  and  in  the  free-stream 

u  =  i)  =  tx;  =  r  =  0,  y  -*  oo;  (19) 

although,  these  can  be  easily  replaced  by  other  conditions  such  as  the  Rankine- 
Hugoniot  conditions  at  the  shock  (Chang  et  al.,  1990)  for  supersonic  flows.  Non- 
homogeneous  boundary  conditions  can  also  be  imposed. 

Nonlinear  PSE 


In  the  linear  PSE  approach  described  above,  the  disturbance  amplitude  is  as¬ 
sumed  to  be  infinitesimally  small  so  that  the  nonlinear  interaction  of  waves  with  dif¬ 
ferent  frequencies  and  spanwise  wave  numbers  is  neglected.  When  finite-amplitude 
waves  are  present  in  the  flow,  the  linear  approach  is  no  longer  valid.  For  nonlinear 
studies,  we  zissume  that  the  total  disturbance  is  again  periodic  in  time  and  in  the 
spanwise  direction,  thus,  the  total  disturbance  function  <j>  can  be  expressed  by  the 
following  Fourier  series 


i(  f*  amn{x)dx+n0z-mu)t) 
*n 


171=— OO  n=— OO 


(20) 


where  a^n  and  'i'mn  are  the  Fourier  components  of  the  streamwise  wave  number 
and  shape  function  corresponding  to  the  Fourier  mode  {mu),n^).  The  frequency 
u  and  wave  number  0  are  chosen  such  that  the  longest  period  and  wave  length 
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are  2:r/u;  and  27r/^  in  the  temporal  and  spanwise  domains,  respectively.  For  most 
stability  problems  of  interest,  it  is  sufficient  to  truncate  (20)  to  only  a  finite  niunber 
of  modes 


M 


N 


<f>=  XI  ^mn(a;,y)e 


«(  f*  am„(i)dz+n^x—mu/t) 
•'*0 


m=  —  M  n=—N 


(21) 


where  M  and  N  zire  the  total  number  of  modes  kept  in  the  truncated  Fourier  series. 
For  all  nonlinear  results  presented  in  this  study,  we  apply  both  the  temporal  and 
spanwise  symmetry  conditions  whenever  applicable,  i.e.,  only  one  quarter  of  modes 
(m  ranging  from  0  to  M  and  n  ranging  from  0  to  N)  are  computed  in  the  marching 
process. 

We  now  substitute  Eq.  (21)  into  our  nonlinear  governing  equation  (6)  and 
perform  harmonic  balance  (collect  terms  with  the  same  spanwise  wave  number  and 
frequency)  for  both  linear  and  nonlinear  terms.  The  resulting  governing  equations 
for  the  shape  function  of  a  single  Fourier  mode  (m,  n)  become 


D  +A  I  n 

■^mn^mn  ^  -^mn  -^mn 


y, 

yy  dy^ 


+  FmnI’A. 


mn 


(22) 


where  matrices  I>mn,  ^mn  and  Bmn  are  given  by 


Dmn  =  +  Z?*  +  iomnA^  +  in^& 

-  -  °^mnWL  + 

Amn  —  A  ^ioCmnV^g  — 

Bjnn  —  B  ioitnnVj.y  —  ifl^Vyg 


and  the  quantity  Amn 


'mn 


»  r *  am„(x)dz 
—  e  *^*0 


The  nonlinear  forcing  function  Fmn  is  the  Fourier  component  of  the  total  forcing 
defined  by  Eq.  (7)  and  can  be  evaluated  by  the  Fourier  series  expansion  of  F" 


Af  N 

F"(x,y,2,t)=  Y,  E  (23) 

m=  —  M  n=  —  N 


The  Fourier  decomposition  of  Eq.  (23)  can  be  done  by  using  the  Fast  Fourier 
Transform  (FFT)  of  F",  which  is  evaluated  numerically  in  the  physical  space.  In 


10 


equation  (22),  a  parabolizing  procedure  similar  to  that  used  in  the  linear  PSE  has 
been  employed  in  order  to  obtain  a  marching  solution. 

As  in  the  linear  PSE,  the  determination  of  the  wave  number  Omn  plays  an 
important  role  in  maintaining  numerical  stability.  The  procedure  described  above 
for  computing  a  for  Unear  disturbances  can  also  be  used  for  the  determination  of 
otmn-  However,  when  all  Fourier  modes  eire  nearly  phase- locked  (as  is  evident  when 
parametric  resonance  of  secondary  instability  takes  place,  see  e.g.  Kachanov  and 
Levchenko  (1984)),  one  may  assume  that  the  wave  number  is  given  as 


^mn  —  ^mn) 


where  Qr  is  the  real  part  of  aio  and  <Tmn  denotes  the  growth  rate  of  the  mode  (m,  n). 
Each  mode  can  have  a  different  imaginary  part  while  the  real  part  is  updated 
according  to  the  phase  change  of  the  dominant  fundamental  mode.  Additional 
phase  shifts  in  the  harmonics  are  included  in  the  evolution  of  the  shape  functions  of 
the  harmonic  waves;  therefore,  all  Fourier  modes  are  not  necessarily  phase-locked. 
It  needs  to  be  pointed  out  that  the  “nearly”  phase-locking  assumption  (since  the 
evolution  of  shape  functions  may  shift  the  phase  slightly)  mentioned  above  is  used 
for  convenience  and  is  not  necessary  for  the  nonlinear  analysis  using  PSE.  In  a 
later  section,  we  will  provide  an  example  of  a  nonUnear  calculation  where  we  let 
the  disturbances  evolve  with  and  without  the  phase-locking  rule.  Use  of  the  phase¬ 
locking  rule,  when  applicable,  saves  computational  cost. 

The  nonlinear  PSE  for  a  single  Fourier  mode,  equation  (22),  is  equivalent  to 
the  linear  PSE  given  in  (14)  with  a  frequency  mu;  and  a  spanwise  wave  number  n/? 
with  the  addition  of  a  forcing  function.  Since  the  forcing  function  acts  eis  a  “source 
term”  of  the  equation,  the  boundary  conditions  and  solution  procedure  described 
above  for  the  linear  PSE  can  be  directly  applied  to  the  nonlinear  system,  except  for 
the  modes  with  zero  frequency  (m  =  0).  These  zero  frequency  modes  are  denoted 
as  the  mean  flow  correction  (if  ra  =  0)  or  longitudinal  vortex  modes  (if  n  ^  0).  For 
these  modes,  as  in  Gortler  vortex  problem,  the  pressure  gradient  dp/dx  drops  out 
making  the  equations  fully  parabolic. 

The  boundary  conditions  given  in  Eqs.  (18)  and  (19)  can  be  applied  to  the 
longitudinal  vortex  mode  without  modification.  For  the  mean  flow  correction,  the 
free-stream  conditions  are  replaced  by 


Woo 


dvoo 

dy 


i&oo  =  Too  =  0) 


y 


oo 


(24) 


to  account  for  the  change  of  displacement  thickness  due  to  the  correction  of  the 
mean  flow  profile  (u  -f  uoo)  arising  from  nonlinear  interactions.  This  Neumann 
condition  for  the  normal  velocity  allows  the  mean  flow  given  by  the  boundary-layer 
solution  to  adjust  itself  in  order  to  assure  mass  balance. 
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Ill  Numerical  Procedure 


In  this  paper,  we  only  consider  the  compressible  stability  of  two-dimensional 
boundary-layer  flow  past  a  flat  plate.  The  mean  flow  solution  is  obtained  by  solv¬ 
ing  the  self-simileir  boundary  layer  equations.  By  using  the  Mangier- Levy- Lees 
transformation,  the  boundary-layer  equations  are  transformed  into  a  set  of  ordi¬ 
nary  differential  equations  (ODE’s).  A  fourth-order  compact  scheme  is  employed  to 
solve  these  ODE’s.  Details  of  the  numerical  procedures  are  given  in  Malik  (1990b) 
and  will  not  be  repeated  here. 

Numerical  solution  of  the  parabolized  stability  equations  (14)  or  (22)  requires 
discretization  in  both  x  and  y  directions.  Since  the  boundary  layer  grows  in  the 
streamwise  direction,  we  expect  that  the  solution  for  the  shape  functions  will  also 
grow.  To  ensure  sufficient  resolution  as  the  disturbances  evolve  downstream,  dis¬ 
cretization  in  the  wall-normal  direction  must  be  able  to  account  for  the  growth  of 
the  boundary  layer.  Instead  of  solving  equations  (14)  and  (22)  in  Cartesian  coor¬ 
dinates,  we  transform  these  equations  to  a  generalized  coordinate  system  defined 
by 

n  =  n{x,y) 

in  order  to  facilitate  numerical  computations  on  a  “growing  mesh”  or  curved  wall 
geometries.  After  this  transformation.  Equation  (14)  becomes 


~d<if  ~d^  - 


where  the  coefficient  matrices  are  given  by 


(26) 


D  =  D 

A  =  ^^A+  ^yB 

V  = 

''VV  —  Vy  ''yy 

The  Jacobian  of  the  transformation  J  is  defined  as 

J  —  ^x^y  ~  ^yVi- 

Equation  (22)  can  be  transformed  in  a  similar  fashion. 

Using  transformation  (25),  we  map  the  computational  grid  into  a  uniform  mesh 
with  constant  increments  in  ^  and  rj  coordinates.  For  most  of  our  calculations,  we 
use  a  constant  step  size  in  x  while  the  grid  is  clustered  near  the  wall  to  resolve  the 
rapid  change  inside  the  boundary  layer.  For  high  Mach  number  calculations,  we  also 
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cluster  the  grid  near  the  critical  layer  located  near  the  boundary-layer  edge.  The 
stretching  along  the  y  direction  is  based  on  the  local  length  scale  /*  (/r  =  y/vexju^) 
which  increases  with  the  boundary-layer  growth.  The  same  grid  distribution  based 
on  yHz  is  used  for  all  x  locations  while  1,  increases  with  x. 

We  use  indices  t  and  j  to  denote  the  grid  index  along  the  streamwise  (x)  and 
wall-normal  (y)  directions,  respectively.  In  the  streamwise  direction,  we  use  the 
second-order  backward  difference 

^ 

for  all  X  locations,  except  for  the  starting  plane  where  a  first-order  backward  differ¬ 
ence  is  employed.  The  resulting  discretized  equation  for  the  i-th  streamwise  plane 
is  then 

~  =  (27) 

In  the  wall-normal  direction,  we  employ  a  fourth-order  accurate  finite-difference 
scheme.  The  two-point  fourth-order  scheme  by  Malik  et  al.  (1982)  is  also  used  for 
normal  derivatives;  however,  it  requires  that  the  normal  mean  velocity  v  is  non-zero 
and  hence  is  not  generally  applicable  for  all  problems. 

For  the  uniform  mesh  in  the  (  —  t}  plane,  the  normal  derivatives  in  Eq.  (26) 
are  discretized  according  to  the  following  fourth-order  central  difference  formulae: 

^  ^  -^i,j+2  -t- 

drj  12At7 

52^ 

^  =(-^i.>+2  + 

-I-  -  ^i,j-2)/l2Arf^. 

For  the  grid  point  next  to  the  boundary,  the  above  five-point  scheme  is  replaced 
with  the  second-order  scheme 


-  2^i,j  +  ^ij-i)/2Ari\ 

At  the  boundary,  five  boundary  conditions  are  needed  for  five  dependent  variables 
in  The  no-slip  and  free-stream  boundary  conditions  given  in  (18)  and  (19) 
are  used.  In  addition,  we  apply  the  discretized  continuity  equation  as  the  fifth 
boundary  condition  both  at  the  wall  and  the  free-stream.  Substituting  the  above 
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normal  derivatives  into  (27)  for  all  interior  points  and  coupling  with  the  boundary 
conditions  result  in  a  block  penta-diagonal  system  of  equations  at  each  x  location 
with  a  block  size  of  5  x  5.  This  block  matrix  can  be  solved  by  the  standard  LU 
decomposition  method. 


IV.  Results  and  Discussion 

To  demonstrate  the  capability  of  the  PSE  approach,  we  perform  both  linear 
and  nonlinesu-  calculations  for  various  Mach  numbers.  In  the  linear  results,  the 
main  focus  will  be  on  the  non-parallel  effect,  and  in  the  nonlinear  regime,  PSE 
calculations  are  carried  all  the  way  to  the  early  stage  of  transition. 

In  the  following  discussion,  we  define  the  growth  rate  in  non-parallel  boimdary 
layers  according  to  Eq.  (17),  i.e.,  for  any  given  flow  variable  ip  (  for  instance,  p,  u, 
etc.),  the  growth  rate  <t  is  defined  as 

a  = -Im(a)  + Re(^^).  (28) 

The  second  term  on  the  right  hand  side  of  the  above  equation  is  a  function  of  y; 
therefore,  the  growth  rate  in  a  non-parallel  boundary  layer  depends  upon  the  dis¬ 
tance  normal  to  the  wall.  We  note  here  that  although  Eqs.  (17)  and  (28)  are  derived 
based  on  the  same  concept,  they  have  different  physical  meanings.  Briefly,  Eq.  (17) 
is  used  to  normalize  the  disturbance  vector  and  determine  the  wave  number  a.  For 
each  normalization,  corresponding  to  different  chosen  in  (17),  the  growth  rate 
for  any  given  variable  at  any  y  location  is  to  be  evaluated  using  Eq.  (28).  For  the 
results  presented  herein,  we  compute  the  growth  rate  at  the  corresponding  loca¬ 
tion  where  the  fluctuation  reaches  its  maximum  value  or  based  on  the  disturbance 
kinetic  energy  integral, 

<TE  =  -/m(o)  +  —{ln\/E) 
ox 

where  E  is  defined  by 

E  =  f  (u^  +v^  +  w^ 

Jo 

for  the  incompressible  limit  and  by 

E  =  I  p{u^  +  +  w^)dy 

Jo 

for  general  compressible  flows.  In  supersonic  wind  tunnel  experiments,  the  growth 
rate  is  usually  measured  for  the  mass  flow  fluctuation.  We  define  the  mass  flow 
fluctuation  as 

(pu)'  =  p'u  -f  pu.  (29) 
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Linear  PSE 


As  mentioned  in  the  previous  section,  the  streamwise  wave  number  a  depends 
upon  the  variable  chosen  and  the  y  location  where  (17)  is  applied.  To  demon¬ 
strate  that  the  resulting  non-parallel  growth  rate  is  very  weakly  dependent  upon 
various  normalizations,  we  first  perform  calculations  for  a  Mach  1.6  boundary  layer 
by  using  different  dependent  variables  to  update  a.  These  variables  include  u,  T 
(evaluated  at  various  y  locations  as  shown  in  the  figure)  and  the  kinetic  energy  inte¬ 
gral  E  defined  above.  Figure  1  shows  the  resulting  imaginary  part  of  the  converged 
value  of  a  for  the  various  norms.  The  results  reveal  that  a  strongly  depends  on 
the  norm  chosen.  The  corresponding  effective  growth  rates  <7^^,  (Tt  and  as,  eval¬ 
uated  by  using  (28)  at  their  maximum  locations  (for  (7p„  and  (Tt  only)  are  shown 
in  Fig.  2.  It  shows  that  the  total  growth  rates  depend  on  how  they  are  measured 
(for  instance,  ar  and  as  are  different);  however,  each  non-parallel  growth  rate  (e.g. 
(Tpu)  appears  to  be  independent  of  the  normalization  procedure  because  results  from 
various  norms  collapse  into  one  single  curve.  Similarly,  adthough  not  shown  here, 
the  non-parallel  wave  number,  evaluated  by 


a 


np 


=  Re(a)  - 


is  also  weakly  dependent  on  the  normalization.  The  above  results  indicate  that 
although  different  norms  result  in  different  values  of  a,  the  total  growth  rate  (and 
wave  number)  by  accounting  for  the  evolution  of  shape  function  in  the  streamwise 
direction  remains  the  same  regardless  of  the  norms.  For  the  results  presented  herein, 
we  use  the  kinetic  energy  integral  E  in  (17)  to  update  a. 

To  verify  the  numerical  algorithm,  the  first  test  case  studied  is  an  incom¬ 
pressible  fiow  ceise.  The  incompressible  results  were  obtained  by  choosing  a  Mach 
number  of  10“®  in  our  compressible  formulation.  Linear  non-parallel  results  are 
available  for  incompressible  boundary  layer  flow  by  using  local  methods  from  many 
authors(e.g.,  Caster,  1974\  The  neutral  points  obtained  from  our  PSE  calcula¬ 
tions  agree  very  well  with  those  from  Caster’s  (1974)  non-parallel  method.  Figure 
3  shows  the  computed  variation  of  the  growth  rates  ((7„,  (t„  and  (Te)  with  Reynolds 
number  {R  =  y/uexfve)  for  a  represen  a'live  non-dimensional  frequency  {F  =  ujR) 
of  1.12  X 10“'* .  The  growth  rates  from  multiple- scales  method  are  shown  by  symbols. 
The  results  shown  in  Figure  3  reveal  that  the  neutral  curve  near  the  upper  branch  is 
shifted  to  higher  Reynolds  numbers  due  to  non-peurallel  effect  as  was  found  by  Caster 
(1974)  and  linear  PSE  results  agree  very  well  with  those  from  the  multiple-scedes 
approach. 

The  second  test  case  was  chosen  to  be  the  Mach  1.6  case  studied  by  El-Hady 
(1991)  using  the  multiple-scales  apn’-'inch.  Tue  frequency  was  fixed  at  0.4  x  10“'* 
and  variable  transport  properties  were  used.  Calculations  were  performed  for  both 
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2-D  and  3-D  linear  disturbzinces  with  an  oblique  wave  angle  of  about  50"  for  the 
latter.  The  growth  rate  of  the  mass  flow  fluctuations  (defined  in  Eqs.  (28)  and 
(29))  from  our  PSE  calculations  together  with  the  multiple- scales  results  are  plotted 
along  with  the  growth  rates  obtained  by  quasi-parallel  linear  stability  theory  in  Fig. 
4.  Our  PSE  results  agree  quite  well  with  those  obtained  from  the  multiple-scales 
approach.  The  results  ailso  indicate  that  for  the  first  mode  disturbance  at  Mach 
1.6,  flow  non-pjirallelism  has  more  effect  on  three-dimensional  disturbances  than 
on  two- dimensioned  ones.  Results  obtained  at  higher  Mach  numbers  also  show  a 
noticeable  non-parallel  effect  on  the  first-mode  instability. 

We  now  show  some  results  for  the  Mach  4.5  flat-plate  flow,  which  is  subject 
to  second-mode  instability  (Mack  1984).  Calculations  were  performed  for  a  dis¬ 
turbance  frequency  of  F  =  1.2  x  10”^  with  different  streamwise  resolutions.  We 
used  step  sizes  ranging  anywhere  from  64  steps  per  wavelength  to  only  one  step 
per  wavelength.  The  results  for  the  second-mode  growth  rate  based  upon  the  to¬ 
tal  kinetic  energy  are  plotted  in  Figure  5.  There  is  essentially  no  difference  in  the 
growth  rate  results  when  two  or  more  steps  per  wavelength  are  used.  The  reason 
why  only  two  points  per  wavelength  could  yield  such  accurate  growth  rates  lies  in 
the  fact  that  most  of  the  wave  information  is  absorbed  in  the  complex  wavenumber 
a.  In  contrast,  direct  numerical  simulation  (DNS)  of  Navier-Stokes  equations  would 
require  many  more  points  per  wavelength  for  comparable  accuracy. 

To  further  verify  our  linear  results,  we  compare  non-parallel  evolution  of  a 
second  mode  disturbance  with  a  frequency  of  2.2  x  10“^  with  DNS.  In  Fig.  6, 
the  maximum  amplitudes  of  various  flow  quantities  are  plotted  agednst  Reynolds 
numbers  for  both  PSE  and  DNS.  The  PSE  results  obtained  by  using  only  7  steps 
per  wave  length  agree  very  well  with  DNS  results  using  16  steps  per  wave  length. 
The  PSE  calculation  took  about  100  seconds  CPU  time  while  the  DNS  required 
more  than  40  hours  on  a  Cray-YMP.  Details  of  the  comparison  including  nonlinear 
disturbances  and  the  spatial  DNS  algorithm  are  given  in  Pruett  and  Chang  (1993). 

Nonlinear  PSE 


a.  Computation  of  a 

In  a  previous  section  we  mentioned  that  the  wavenumber  a  for  the  harmonics 
may  be  determined  either  by  the  phase-locking  rule  or  by  using  Eq.  (17).  It  is  known 
that  the  nonlinear  wave  interaction  is  dependent  on  the  phase- difference  between 
various  modes.  Therefore,  it  is  essential  that  nonlinear  PSE  approach  must  not 
require  phase-locking  as  a  fundamental  assumption;  although,  it  may  be  used  as  a 
convenience  for  problems  where  phase- locking  happens  anyway. 

In  order  to  demonstrate  that  phase-locking  is  not  a  basic  assumption  for  non¬ 
linear  PSE  computations,  the  following  test  has  been  performed.  Nonlinear  calcula¬ 
tions  have  been  done  for  a  flat-plate  boundary  layer  in  the  incompressible  limit.  A 
two-dimensional  wave  with  a  frequency  F  =  .86  x  10~^  and  an  initial  amplitude  of 
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.25%  at  /2  =  400  is  introduced  in  the  boundary  layer  and  the  evolution  of  this  wave 
along  with  its  various  harmonics  is  monitored.  Calculations  were  performed  in  two 
different  ways.  First,  the  wavenumber  a  was  computed  for  the  fundamental  wave 
according  to  Eq.  (17)  and  the  phase-locking  rule  was  used  for  all  the  harmonics.  In 
the  second  set  of  calculations,  wavenumbers  for  the  fundamental  and  all  the  harmon¬ 
ics  were  computed  independently  by  using  Eq.  (17).  The  computed  results  for  the 
amplitude  of  u  velocity  (fundamental,  its  four  harmonics  and  meanflow  correction) 
are  presented  in  Fig.  7(a).  It  can  be  seen  that  only  very  minor  differences  appear 
between  the  two  sets  of  calculations  and  these  also  tend  to  disappear  as  the  calcu¬ 
lations  are  marched  away  from  the  inflow  boundary.  The  second  set  of  calculations 
takes  about  50%  more  computer  time  due  to  the  iterations  involved  in  determining 
a  for  the  wave  harmonics.  Hence,  it  is  expedient  to  use  the  phase-locking  rule  for 
problems  where  this  may  be  the  outcome  in  any  case. 

In  Fig.  7(b),  the  same  nonlinear  PSE  results  are  compared  with  the  spatial 
incompressible  DNS  results  for  fundamental,  first  harmonic  and  the  mean  flow  dis¬ 
tortion  modes.  Both  PSE  and  DNS  start  with  the  same  initial  conditions,  i.e.,  a  fun¬ 
damental  disturbance  (1,0)  at  R  =  400  and  all  harmonic  waves  including  the  mean 
flow  distortion  are  generated  through  nonlinear  interactions.  The  good  agreement 
between  DNS  and  nonlinear  PSE  indicates  that  the  parabolizing  approximation  in 
the  PSE  approax:h  does  not  introduce  any  severe  error  and  all  detailed  nonlinear 
features  are  properly  captured.  Details  of  the  comparison  including  disturbance 
profiles  can  be  found  in  Joslin  et  al.  (1992). 
b.  Second-Mode  Instability  at  Mach  4.5 

To  verify  the  nonlinear  PSE  algorithm,  we  choose  the  nonlineeu:  second  mode 
simulation  at  Mach  4.5  investigated  by  Erlebacher  and  Hussaini  (1990)  using  the 
temporal  DNS  approach.  As  in  the  temporal  DNS  approach,  we  assume  that  the 
mean  flow  is  partdlel  and  study  the  spatial  evolution  of  disturbances  in  the  presence 
of  nonlinear  interactions.  The  initial  conditions  were  provided  by  the  eigensolution 
from  the  linear  theory  at  J?  =  781  and  four  Fourier  modes  (M  =  3)  were  kept  in  the 
truncated  series.  In  our  PSE  calculation,  the  disturbance  is  assumed  to  be  periodic 
in  time  and  the  nonlinear  evolution  is  carried  downstream  in  z  as  opposed  to  the 
temporal  DNS  approach  where  the  disturbance  is  periodic  in  z  and  integration  is 
carried  in  time.  It  was  found  in  Erlebacher  and  Hussaini  (1990)  that  due  to  non¬ 
linear  effect,  the  growth  rate  of  the  fundamental  disturbance  strongly  depends  on 
y  and  there  exists  a  sharp  decrease  in  the  local  growth  rate  near  the  critical  layer. 
The  growth  rates  based  on  «io  from  our  PSE  results  are  shown  in  Figure  8(a)  for 
different  z  locations.  The  growth  rate  is  initially  uniform  at  the  starting  location 
(z  =  OA).  As  the  disturbances  are  evolving  downstreeim,  nonlinear  effects  observed 
in  Erlebacher  and  Hussaini  (1990)  are  evident  in  the  present  spatial  calculations. 
For  comparison,  their  temporal  DNS  results  are  shown  in  Figure  8(b)  at  different 
time  levels  represented  as  multiples  of  the  temporal  period  r.  Figure  9  depicts  the 
amplitudes  of  the  density  fluctuation  of  the  first  harmonic  for  both  PSE  calcula- 


17 


tions  and  the  DNS  results.  The  DNS  results  in  Figure  9  are  re-scaled  to  facilitate 
comparison.  Qualitatively,  all  nonlinear  features  observed  in  DNS,  including  the 
kink  near  the  boundary-layer  edge,  are  properly  resolved  in  our  PSE  results, 
c.  Subharmonic  and  Fundamental  Resonance  at  Mach  1.6 

Numerical  simulation  of  incompressible  flows  have  shown  that  the  rapid  growth 
of  three-dimensional  secondary  disturbances  is  followed  by  breakdown  to  turbulence. 
Secondary  instability  is  triggered  when  the  primary  disturbances  reach  sufficiently 
high  amplitudes.  To  show  the  capability  of  the  PSE  approach  in  simulating  transi¬ 
tion  onset,  we  also  perform  nonlinear  calculations  to  study  the  secondary  instability 
mechanism.  We  carry  our  calculations  all  the  way  to  transition  for  both  K-type  (fun¬ 
damental)  and  H-type  (subharmonic)  breakdown.  We  choose  Mach  1.6  flat  plate 
flow  with  a  primary  disturbance  frequency  of  0.5  x  10“'*.  The  same  flow  conditions 
were  also  used  by  Thumm  et  al.  (1989)  in  their  spatial  Navier-Stokes  simulations 
of  a  compressible  boundary  layer. 

We  first  perform  a  series  of  calculations  to  determine  the  amplitude  of  the 
primary  disturbance  which  will  trigger  the  secondary  instability.  The  PSE  calcu¬ 
lation  is  initiated  at  a  Reynolds  number  R  of  460  where  we  impose  a  primary  2-D 
wave  (mode  (2, 0))  obtained  by  a  local  eigenvalue  calculation  and  two  subharmonic 
waves  ((1,1)  and  (1,-1)  modes)  by  using  the  compressible  secondary  instability 
theory  (Ng  and  Erlebacher,  1992)  and  all  the  remaining  harmonics  are  assumed  to 
have  zero  amplitudes.  Initial  amplitudes  of  the  primary  disturbances  are  set  to  be 
3%,  1.1%  and  0.6%  at  the  inflow  plane  which  corresponds  to  the  5%,  2%  and  1% 
(the  maximum  amplitudes  near  the  vibrating  ribbon)  cases  given  in  Thumm  et  al. 
(1989).  The  initial  amplitudes  of  the  subharmonic  waves  are  fixed  at  0.019%  for 
all  three  cases.  The  spanwise  wave  number  of  the  subharmonic  mode  is  fixed  at 
^/R  =  0.53  X  10”*  which  corresponds  to  an  oblique  wave  angle  of  45°.  Six  temporal 
Fourier  modes  and  three  spanwise  modes  (M  =  5  and  N  =  2)  are  kept  in  the  Fourier 
series.  The  evolution  of  both  primary  and  subharmonic  disturbances  are  shown  in 
Figure  10.  Qualitatively,  our  results  agree  with  those  of  Thumm  et  al.(1989).  Any 
quantitative  differences  are  due  to  different  initial  conditions.  We  find  that  a  1.1% 
initial  amplitude  (2%  case  in  Thumm  et  al.  (1989))  for  the  primary  mode  is  enough 
to  trigger  the  secondary  growth;  however,  the  onset  of  secondary  growth  for  this 
case  occurs  at  i?  =  800  where  the  primary  wave  is  about  to  decay.  We  continue 
the  PSE  calculations  beyond  R  =  1050  for  this  case  and  find  that  the  secondary 
disturbance  eventually  saturates  and  the  flow  does  not  reach  the  transitional  stage. 

To  carry  the  3%  case  to  the  transition  stage,  we  made  another  calculation  with 
more  Fourier  modes  (M  =  7  and  iV  =  4)  and  the  maximum  velocity  amplitudes  of 
some  representative  modes  aue  given  in  Figure  11.  Besides  the  fundamental  mode 
(2,0)  and  the  subharmonic  mode  (1,1),  higher  harmonics  in  time  aind  spanwise 
domain  are  also  excited  due  to  nonlineau*  interaction.  Initially,  the  (4, 0)  mode  gains 
energy  from  self-interaction  of  the  (2, 0)  mode  and  the  interaction  of  (2, 0)  aind  (1,1) 
produces  the  (3, 1)  mode.  When  the  subhairmonic  mode  grows  due  to  the  onset  of 
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secondary  instability,  its  harmonic  (2,2)  also  grows  at  slightly  higher  rate.  The 
streamwise  vortex  mode  (0,2)  arises  due  to  the  interaction  of  (1,1)  mode  zuid  its 
complex  conjugate  (—1,1).  As  all  these  modes  continue  to  grow,  more  eind  more 
modes  are  excited.  The  energy  cascade  exhibits  a  staggered  pattern.  For  instance, 
among  the  two-dimensional  modes,  only  (2,0),  (4,0),  (6,0),  etc.  gain  energy;  while 
for  modes,  only  (1,1),  (3, 1),  (5, 1),  etc.  are  excited.  The  remaining  modes  (e.g. 
(1,0),  (3,0),  (0,1),  (2,1)  etc.)  remain  unexcited  throughout  the  calculation.  The 
above  staggered  energy  cascade  is  typical  for  subharmonic  secondary  instability. 
The  secondary  amplitude  overtakes  the  primary  at  about  R  =  980  and  reaches  an 
equilibrium  state  around  R  =  1100.  At  this  stage,  many  harmonic  waves  reach 
fairly  high  amplitudes  as  the  flow  heads  for  transition.  We  plot  the  time  sequence 
of  spanwise  vorticity  contours  at  the  peak  and  valley  planes  (corresponding  to  the 
maximum  and  minimum  disturbance  rms  amplitudes,  respectively)  in  Figures  12(a) 
and  12(b).  As  can  be  seen,  the  vorticity  pattern  doubles  its  wavelength  for  i  > 
2200  (a:  is  normalized  w.r.t.  the  boundary  layer  length  scale  I  at  the  initial  plane) 
indicating  the  presence  of  high- amplitude  subhaxmonic  wave.  It  is  also  evident 
that  the  vortex  roll-up  results  in  a  distinct  kink  in  the  shear  layer.  Towards  the 
end  of  the  computation,  regions  of  intense  vorticity  near  the  wall  begin  to  appear 
indicating  that  flow  is  heading  for  breakdown.  Figure  13  shows  the  streamwise 
velocity  contours  in  the  x-z  plane  for  a  wall  normal  distance  of  y  =  2.3,  where  the 
TS  wave  reaches  its  maximum  according  to  the  linear  solution.  The  flow  is  initially 
two-dimensional  and  three-dimensional  effect  becomes  important  for  x  >  1800.  For 
X  >  2000,  a  staggered  contour  pattern  is  evident.  This  pattern  is  associated  with 
the  lambda  vortex  structure,  a  distinct  characteristic  of  subharmonic  breakdown, 
as  observed  in  many  incompressible  experiments  (e.g.  Corke  and  M^mgano  (1989)). 

Nonlinear  PSE  calculations  are  also  performed  for  the  same  Mach  1.6  case  but 
for  a  fundamental- type  secondary  resonance.  The  initial  amplitude  of  the  primary 
wave  is  again  3%  and  that  of  the  secondary  is  taken  to  be  0.005%  to  minimize 
nonlinear  interaction  close  to  the  starting  location.  The  spanwise  wave  number  is 
j3fR  =  1.52  X  10“"*  (oblique  wave  angle  of  60®  for  the  secondary  wave)  and  the  pri¬ 
mary  wave  frequency  is  again  0.5  x  10““*.  The  initial  conditions  for  our  marching  cal¬ 
culation  consist  of  a  2-D  primary  wave  (mode  (1, 0)),  two  oblique  fundamental-type 
secondary  disturbances  (mode  (1,1),  (1,-1))  and  the  longitudinal  vortex  (mode 
(0, 1)).  The  same  number  of  Fourier  modes  as  in  the  subharmonic  case  is  used. 

Nonlinear  evolution  of  the  maximum  rms  amplitude  of  u'  is  shown  in  Figure  14. 
Initially,  the  dominant  modes  eire  (1, 0),  (1, 1),  (0, 1)  and  (2, 0)  (the  first  harmonic  of 
the  fundamental  2D  mode).  Unlike  the  subharmonic  case,  all  harmonic  waves  (both 
odd  and  even  modes)  gain  energy  directly  from  nonlinear  interaction.  Among  them, 
the  (2,1)  (due  to  (1,0)  and  (1,1))  and  (1,2)  (due  to  (0,1)  and  (1,1))  modes  are 
more  noticeable.  For  Reynolds  numbers  beyond  870,  the  spectrum  is  rapidly  filled 
with  high-amplitude  disturbances  and  the  flow  is  heading  towards  transition.  As 
compared  to  the  subharmonic  breakdown,  transition  location  shifts  upstream  due 
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to  the  larger  growth  rate  of  the  secondary  disturbance  as  a  consequence  of  higher 
oblique  wave  angle. 

The  time  sequence  of  spanwise  vorticity  contours  over  a  period  of  the  primary 
wave  is  shown  in  Figures  15(a)  and  15(b)  for  the  peak  and  valley  planes,  respec¬ 
tively.  In  contrast  to  the  subharmonic  breakdoown,  the  wave  length  remains  the 
same  throughout  the  whole  computational  dommn.  One  important  characteristic  of 
the  K-type  breakdown  is  the  appearance  of  aligned  lambda  vortices.  This  is  better 
visualized  in  the  streamwise  velocity  contours  shown  in  Figure  16  for  x  >  1300. 
Similar  to  that  observed  in  incompressible  simulations  of  Zsing  and  Krist  (1989), 
regions  of  intense  shear  begin  to  appear  near  the  end  of  the  computational  dommn 
in  Figures  15(a)  and  15(b).  This  indicates  that  flow  has  just  entered  the  transitional 
stage.  It  is  confirmed  by  plotting  the  average  wall  shear  in  Figure  17.  The  com¬ 
puted  wall  shear  is  slightly  above  the  laminar  value  for  most  of  the  computational 
domain.  Only  towards  the  end,  wall  shear  significantly  departs  from  the  laminm 
value  indicating  the  onset  of  transition.  In  this  way,  PSE  provides  the  prediction 
of  boundary-layer  transition  for  the  imposed  initial  conditions.  The  PSE  wjdl  shear 
lies  above  the  laminar  value  right  from  the  beginning  because  of  the  relatively  high 
amplitude  t«f  the  2-D  primary  disturbance  needed  for  transition  in  supersonic  fiow. 
Since  most  amplified  waves  in  supersonic  flow  are  not  two-dimensional,  oblique 
primary  modes  may  lead  to  transition  for  lower  initial  amplitudes.  In  order  to 
carry  the  calculations  further  into  transitional  regime,  more  spanwise  and  temporal 
resolution  will  be  required.  It  remains  to  be  seen  how  far  PSE  can  proceed  into 
the  transitional  zone.  The  computational  time  used  for  the  results  presented  in 
Figures  14-17  was  15  minutes  on  a  Cray-YMP  machine.  Similar  results  from  full 
compressible  Navier-Stokes  equations  would  require  0(50)  hours. 

V.  Conclusions 

Linear  and  nonlinear  compressible  boundary-layer  stability  is  studied  by  using 
the  PSE  approach.  Several  issues  concerning  the  characteristics  of  the  paraboliza¬ 
tion  and  the  updating  of  the  streamwise  wave  number  are  also  discussed.  The 
governing  equations  are  solved  by  using  second-order  backward  differences  for  the 
streamwise  derivatives  while  the  wall-normal  direction  is  discretized  by  a  fourth- 
order  accurate  finite-difference  scheme. 

Non-parallel  fiow  effects  have  been  studied  for  linear  disturbances.  For  oblique 
waves  of  the  first  mode  type,  the  departure  from  the  parallel  results  is  more  pro¬ 
nounced  as  compared  to  that  for  the  two-dimensional  waves.  Our  linear  results  are 
in  good  agreement  with  those  from  the  multiple-scales  approach,  as  well  as  those 
from  full  Navier-Stokes  equations. 

Nonlinear  PSE  calculations  are  carried  all  the  way  to  the  early  stage  of  tran¬ 
sition  for  a  Mach  1.6  flow.  Both  the  subharmonic  and  fundamental  types  of  break- 
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down  axe  studied  by  the  current  PSE  approach.  Qualitatively,  these  breakdown 
processes  are  similar  to  the  ones  in  incompressible  boundary  layers,  except  that 
high  amplitudes  of  the  2-D  primary  wave  are  required.  The  promising  results  of 
our  PSE  calculations  show  that  this  new  approach  is  a  powerful  tool  for  the  study 
of  boundary-layer  stability  and  transition  prediction.  The  parabolized  form  of  the 
governing  equations  allow  the  numerical  solution  to  be  obtained  in  a  computational 
time  which  is  orders  of  magnitude  lower  than  that  required  for  direct  simulation  of 
Navier-Stokes  equations. 


References 

Bertolotti,  F.  P.,  Herbert,  Th.  and  Sp2dart,  P.  R.  1992  “Linear  and  Nonlinear 
stability  of  the  Blasius  Boundary  Layer,”  J.  Fluid  Mech.,  Vol.  242,  pp.441-474. 

Chang,  C.-L.,  Malik,  M.  R.,  and  Hussaini,  M.  Y.  1990  “Effects  of  Shock  on  the 
Stability  of  Hypersonic  Boundary  Layers,”  AIAA  Paper  No.90-1448. 

Chang,  C.-L.  and  Merkle,  C.  L.  1989  “The  Relation  between  Flux  Vector  Split¬ 
ting  and  Parabolized  Schemes,”  J.  Computational  Physics,  Vol.  80,  No.  2,  pp. 
344-361. 

Corke,  T.  C.  and  Mangano,  R.  A.  1989  “Resonant  Growth  of  Three-dimensional 
Modes  in  Transitioning  Blausius  Boundary  Layers,”  J.  Fluid  Mech.,  Vol.  209,  pp.93- 
150. 

Denier,  J.  P.,  Hall,  P.,  and  Seddougui,  S.  O.  1991  “On  the  Receptivity  Problem 
for  Gortler  Vortices  and  Vortex  Motions  Induced  by  Wall  Roughness,”  Phil.  Trans. 
R.  Soc.  Lond.  A,  335,  pp. 51-85. 

El-Hady,  N.  M.  1991  “Nonparallel  Instability  of  Supersonic  and  Hypersonic 
Boundary  Layers,”  Phys.  Fluids  A,  Vol.  3,  No.9,  pp.  2164-2178. 

Erlebacher,  G.  and  Hussaini,  M.  Y.  1990  “Numerical  Experiments  in  Supersonic 
Boundary- Layer  Stability,”  Phys.  Fluids  A,  Vol.2  No.l,  pp.94-104. 

Gaponov,  S.  A.  1981  “The  influence  of  Flow  Non-Parallelism  on  Disturbance 
Development  in  the  Supersonic  Boundary  Layers,”  Proceedings  of  the  8-th  Canadian 
Congress  of  Applied  Mechanics,  pp.  673-674. 

Gaster,  M.  1974  “On  the  Effects  of  Boundary-Layer  Growth  on  Flow  Stability,” 
J.  Fluid  Mech.,  Vol.  66,  pp.465-480. 

Hall,  P.  1983  “The  Linear  Development  of  Gortler  Vortices  in  Growing  Bound¬ 
ary  Layers,”  J.  Fluid  Mech.,  Vol.  130,  pp.41-58. 

Joslin,  R.  D.,  Streett.  C.  L.  and  Chang,  C.-L.  1992  “Validation  of  Three- 
Dimensional  Incompressible  Spatial  Direct  Numerical  Simulation  Code  -  A  Compar¬ 
ison  with  Linear  Stability  and  Parabolic  Stability  Equation  Theories  for  Boundary- 
Layer  Transition  on  a  Flat  Plate,”  NASA  Technical  Paper  TP-3205. 

Kachanov,  Y.  S.  and  Levchenko,  V.  Y.  1984  “The  Resonant  interaction  of 
disturbances  at  Laminar- turbulent  transition  in  a  Boundary  Layer,”  J.  Fluid  Mech., 


21 


Vol.  138,  pp.209-247. 

Mack,  L.  M.  1969  “Boundary-layer  Stability  Theory,”  Document  900-277,  Rev. 
A,  JPL,  Pasadena. 

Mack,  L.  M.  1984  “Boundary- Layer  Linear  Stability  Theory,”  AGARD  Report 
No.  709. 

Malik,  M.  R.,  Chuang,  S.,  and  Hussaini,  M.  Y.  1982  “  Accurate  Numerical 
Solution  of  Compressible,  Linear  Stability  Equations,”  ZAMP,  Vol.  33,  189. 

Malik,  M.  R.  1990a  “Stability  Theory  for  Laminar  Flow  Control  Design,”  Vw- 
cous  Drag  Reduction  in  Boundary  Layers,  ed.  D.  M.  Bushnell  and  J.  N.  Hefner, 
pp.3-46. 

Malik,  M.  R.  1990b  “Numerical  Methods  for  Hypersonic  Boundary  Layer  Sta¬ 
bility,”  J.  Computational  Physics,  Vol.  86,  No.  2,  pp. 376-413. 

Ng,  L.  and  Erlebacher,  G.  1992  “Secondary  Instabilities  in  Compressible 
Boundary  Layers”,  Phys.  Fluids  A,  Vol.  4,  No.  4,  pp. 710-726. 

Pruett,  C.  D.  and  Chang,  C.-L.  1993  “A  Comparison  of  PSE  and  DNS  for  High- 
Speed  Boundary-Layer  Flows,”  ASME  Fluids  Engineering  Conference,  Washington, 
D.  C.,  June  21-24,  1993. 

Rubin,  S.  G.  1981  “A  Review  of  Marching  Procedures  for  Parabolized  Navier- 
Stokes  Equations,”  Proceedings  of  Symposium  on  Numerical  and  Physical  Aspects 
of  Aerodynamic  Flows,  Springer- Verlag,  New  York,  pp. 171-186. 

Saric,  W.  S.  and  Nayfeh,  A.  H.  1975  “Non-parallel  Stability  of  Boundary- Layer 
Flows,”  Phys.  Fluids,  Vol.  18. 

Spall,  R.  E.  and  Malik,  M.  R.  1989  “Gortler  Vortices  in  Supersonic  zmd  Hy¬ 
personic  Boundary  Layers,”  Phys.  Fluids  A,  Vol.  1,  No.ll,  pp.  1822-1835. 

Thumm,  A,  Wolz,  W.,  and  Fasel,  H.  1989  “Numerical  Simulation  of  Spatially 
Growing  Three-Dimensional  Disturbance  Waves  in  Compressible  Boundary  Lay¬ 
ers,”  Proceedings  of  the  third  I U TAM  Symposium  on  Laminar- Turbulent  Transition, 
Toulouse,  France,  Sept.  11-15,  1989. 

Vigneron,  Y.  C.,  Rakich,  J.  V.,  and  Tannehill,  J.  C.  1978  “Calculation  of 
Supersonic  Viscous  Flow  over  Delta  Wings  with  Sharp  Subsonic  Leading  Edges,” 
AIAA  Paper  No.  78-1337. 

Zang,  T.  A.  and  Krist,  S.  E.  1989  “Numerical  Experiments  on  Stability  and 
Transition  in  Plane  Channel  Flow,”  Theoretical  &  Comput.  Fluid  Dynamics,  Vol. 
1  No.  1,  pp. 41-64. 


22 


0.00200 


"  /y 

^  y  ^ 

/  y  y 

■'  /  y 
y  /  y 

'  /  //  / 

/  ///  / 

i///' 

/  r//  / 


/  # 

i/l 


11*  CO  ll*  CVJ  CO 

>»  JL  4L  iL  "  “ 

ULJ  jjS  .4^  <4^  - 

C3  05  C3  Cfl  ^ 

C  3  ^  3  4-. 


V^\^ 

\\\ 

\\\ 

\n 


\\  V. 

\  \  \  N 

\  \  \  \ 
V  \  • 

\\ 


Figure  1.  converged  values  of  —Oi  by  using  various  normalizations  for  a  Mach  1.6 
flow  with  F  =  0.5  X  10“^ 
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incompressible  nonlinear  F=86 
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Figure  7.  Nonlinear  evolution  of  a  fundamental  2-D  wave  with  F  =  0.86  x  10“^  in 
the  incompressible  limit:  (a)8howing  the  velocity  amplitudes  both  with  and  without 
phase-locking  rule  (b)comparing  with  spatial  DNS  results. 
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Mach  4.5  Non-linear  Second  Mode  R*781 
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Figure  8.  Nonlinear  second  mode  simulation  on  a  parallel  mean  flow  of  Mach  4.5 
and  i?  =  781,  showing  the  evolution  of  growth  rates  against  y  at  various  downstream 
distances  represented  as  multiples  of  the  fundamental  wave  length  A  (for  PSE)  or 
the  temporal  period  r  (for  DNS).  (a)Nonlinear  PSE  (b)DNS 
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Mach  4.5  Non-linear  Second  Mode  R=781 


Figure  9.  Nonlinear  second  mode  simulation  on  a  parallel  mean  flow  of  Mach  4.5  and 
R  —  781:  showing  the  evolution  of  density  profiles  against  y  for  the  first  harmonic 
wave  (mode  (2,0)  ).  Lines  are  from  nonlinear  PSE  calculations  and  symbols  are 
from  DNS. 


Subharmonic  Breakdown  F=50 


Figure  11.  Amplitude  evolution  of  various  Fourier  modes  for  a  Mach  1.6  subhar¬ 
monic  breakdown.  (Initial  amplitudes  for  the  primary  and  subharmonic  modes  are 
the  same  as  the  3%  case  in  Figure  8). 


Figure  12.  Time  sequence  of  spanwise  vorticity  contours  for  the  Mach  1.6  subhar¬ 
monic  breakdown  (x  and  y  are  streamwise  and  wall-normal  coordinates  normalized 
by  the  boundary-layer  length  scale  at  i?  =  460):  (a)  Peak  plane,  (b)  Valley  plane. 
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Subharmonic  Breakdown 
Dity  Contours  at  y/lo=2.3 


ach  1.6  Fundamental  Breakdown  F=50 


Figure  14.  Amplitude  evolution  of  various  Fourier  modes  for  a  Mach  1.6  fundamen¬ 
tal  breakdown  (F  =  0.5  x  10~^  and  ^/R  =  1.52  x  10“^)  with  initial  amplitudes  of 
3%  and  0.005%  for  the  primary  and  secondary  disturbances  at  il  =  460. 


Figure  15.  Time  sequence  of  spanwise  vorticity  contours  for  the  Mach  1.6  funda¬ 
mental  breakdown  :  (a)  Peak  plane,  (b)  Valley  plane. 
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Fundamental  Breakdown 
jity  Contours  at  y/l^=2.3 
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Figure  16.  Strecunwise  velocity  contours  in  x—z  plane  at  y  =  2.3  for  the  fundamental 
breakdown. 
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